model {
  for (i in 1:NOBS) {
    #Stochastic component
    s[i,1:NPARTY] ~ dmulti(pi[i,1:NPARTY], n[i])

    #Systematic component
    for (j in 1:NPARTY) {
      mu[i,j] <- lambda[j] + rho*log(v[i,j])
      expmu[i,j] <- exp(mu[i,j])
      pi[i,j] <- expmu[i,j]/sum(expmu[i,1:NPARTY])
    }
  }

  ## Priors
  rho ~ dunif(2,3) #Uniform
  lambda[1] <- 0 # Restriction
  for (j in 2:NPARTY) {
    lambda[j] ~ dnorm(0, 1.0E-4) # Flat normal priors
  }

  ## Expected an Predicted values
  for (j in 1:NPARTY) {
    xb[j] <- lambda[j] + rho*log(xstar)
    exb[j] <- exp(xb[j])
    ev[j] <- exb[j]/sum(exb[1:NPARTY]) #Expeced values (Probabilities)
  }
  pv[1:NPARTY] ~ dmulti(ev[1:NPARTY], 300) #Predicted values

  ## Model fit
  for (i in 1:NOBS) {
    Pred[i,1:NPARTY] ~ dmulti(pi[i,1:NPARTY], n[i])
  }
}